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Viscosity  parameterization  and  the  Gulf  Stream 
separation 

Eric  P.  Chassignet  and  Zulema  D.  Garraffo 

Rosenstiel  School  of  Marine  and  Atmospheric  Science,  University  of  Miami,  Miami,  Florida  - USA 


Abstract.  Recent  advances  in  computer  architecture  allow  for  numerical  in- 
tegration of  state-of-the-art  ocean  models  at  basin  scale  with  a grid  resolution 
of  1/10°  or  higher.  At  that  resolution,  the  Gulf  Stream’s  separation  at  Cape 
Hatteras  is  well  simulated,  but  substantial  differences  from  observations  are 
still  observed  in  its  path,  strength,  and  variability.  Several  high  resolution 
(1/12°)  North  Atlantic  simulations  performed  with  the  Miami  Isopycnic  Co- 
ordinate Ocean  Model  (MICOM)  are  discussed  and  the  results  suggest  that, 
even  with  such  a fine  grid  spacing,  the  modeled  large  scale  circulation  is  still 
quite  sensitive  to  choices  in  forcing  and  viscosity  parameterization. 


1.  Introduction 

Until  recently,  most  ocean  general  circulation  models 
(OGCMs)  had  great  difficulties  in  reproducing  the  basic 
pattern  of  the  Gulf  Stream.  The  modeled  Gulf  Stream 
had  in  general  the  tendency  to  separate  far  north  of 
Cape  Hatteras  and  to  form  a large  stationary  anticy- 
clonic  eddy  at  the  separation  latitude  [see  Dengg  et  al. 
(1996)  for  a review].  Simulations  with  grid  resolution  of 
1/10°  or  higher  are  now  able  to  realistically  represent 
the  Gulf  Stream  separation  ( Paiva  et  al .,  1999;  Smith 
et  al.,  2000;  Hurlburt  and  Hogan,  2000).  These  results 
support  the  view  that  a good  representation  of  the  in- 
ertial boundary  layer  is  an  important  factor  in  the  sep- 
aration process  ( Ozgdkmen  et  al.,  1997).  The  fine  mesh 
size  also  resolves  the  first  Rossby  radius  of  deforma- 
tion everywhere  in  the  subtropical  gyre  (marginally  in 
the  subpolar  gyre),  therefore  providing  a good  represen- 
tation of  baroclinic  instability  processes  ( Paiva  et  al., 
1999;  Smith  et  al.,  2000). 

However,  despite  the  more  realistic  behavior,  the  rep- 
resentation of  the  Gulf  Stream  separation  differs  from 
one  simulation  to  the  next,  sometimes  significantly. 
This  paper  discusses  some  of  the  factors  influencing  the 
modeled  circulation  in  the  Miami  Isopycnic  Coordinate 
Ocean  Model  (MICOM).  It  will  be  shown  that  even  with 
such  a fine  grid  spacing,  the  viscosity  parameterization 
remains  of  importance  for  the  modeled  large  scale  ocean 
circulation. 

2.  Mean  sea  surface  height  fields 

Figures  1 and  2 display  the  mean  sea  surface  height 
(SSH)  from  two  MICOM  simulations  with  an  horizon- 


Figure  1.  3-year-mean  model  SSH  field  for  the  1/12° 
COADS-forced  MICOM.  The  viscosity  operator  is  a 
combination  of  biharmonic  ( At  = VpAx3,  with  Vfi  = 1 
cm/s)  and  Laplacian  (A i = max  [.1  Ax2x  deformation 
tensor,  VpAar],  with  Vp  = .5  cm/s). 


tail  resolution  of  A <f>  = 1/12°.  The  horizontal  grid  is 
defined  on  a Mercator  projection  with  the  resolution 
given  by  A <j>  x A^cos(^),  where  <j>  is  the  latitude.  The 
first  simulation  (Fig.  1),  configured  from  28°S  to  65°N, 
was  integrated  with  MICOM  for  20  years  using  monthly 
climatological  COADS-based  forcing  (including  fresh- 
water flux)  plus  a weak  restoration  to  monthly  climato- 
logical surface  salinity  ( Paiva  et  al.,  1999,  Garraffo  et 
al.,  2001a, b).  The  second  simulation  (Fig.  2),  config- 
ured from  28°S  to  70°N,  including  the  Mediterannean 
Sea,  was  first  spun-up  for  6 years  with  MICOM  us- 
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Figure  2.  2-year-mean  model  SSH  field  for  the  1/12° 
ECMWF-forced  MICOM.  The  viscosity  operator  is  the 
same  as  for  the  COADS-forced  run  (see  caption  of 
Fig.  1). 

ing  monthly  climatological  ECMWF  atmospheric  fields 
(including  freshwater  flux)  plus  a weak  restoration  to 
monthly  climatological  surface  salinity,  and  is  presently 
further  integrated  using  6-hourly  ECMWF  forcing  from 
1979  to  2000. 

In  both  simulations,  the  simulated  Gulf  Stream  path 
agrees  well  with  observations  until  the  location  of  the 
New  England  Seamounts  chain.  Eastward  of  the  chain, 
the  ECMWF-forced  run  (Fig.  2)  exhibits  a path  that 
agree  well  with  observations  everywhere.  In  the  COADS- 
forced  run  (Fig.  1),  the  path  east  of  the  New  England 
Seamounts  chain  is  displaced  northward  by  about  1°  to 
2°.  This  northward  shift  in  the  COADS-forced  run  is 
associated  with  a larger  than  observed  seasonal  migra- 
tion of  the  path  [observed  annual  signal  of  up  to  100 
km,  north  of  the  mean  from  August  to  November  and 
south  of  the  mean  from  March  to  June  (A.  Mariano, 
1999,  personal  communication)].  This  higher  than  ob- 
served seasonal  shift  results  primarily  from  the  fact  that 
the  MICOM  bulk  Kraus-Turner  mixed  layer  is  on  the 
average  deeper  in  the  COADS-forced  run  than  in  the 
ECMWF-forced  run  (not  illustrated).  The  deepening  of 
the  mixed  layer  in  winter  induces  a decrease  in  the  mag- 
nitude of  the  upper  layer  velocities  because  MICOM’s 
mixed  layer  does  not  allow  vertical  shear  (Kraus  and 
Turner,  1967).  A deeper  mixed  layer  in  the  COADS- 
forced  run  therefore  implies  a Gulf  Stream  that  is  less 
inertial  than  in  the  ECMWF-forced  run.  The  end  re- 
sult is  that  in  the  latter  run,  the  modeled  Gulf  Stream 
path  agrees  well  with  observations  and  does  not  exhibit 
a higher  than  observed  seasonal  shift  in  latitude  east- 
ward of  the  New  England  seamounts  chain. 


The  impact  of  the  seamounts  on  the  Gulf  Stream 
path  and  variability  was  further  investigated  in  a 3-year 
sensitivity  experiment  with  COADS  forcing  in  which 
the  bottom  topography  was  modified  by  removing  the 
New  England  seamounts.  The  impact  of  removing  the 
seamounts  on  the  Gulf  Stream  path  was  found  to  be 
negligible  (not  illustrated). 

3.  Importance  of  the  viscosity 
parameterization 

When  the  grid  spacing  reaches  a certain  threshold, 
the  energy  cascade  from  the  small  to  the  large  scales 
should  be  properly  represented  by  the  model  physics. 
Dissipation  should  then  be  prescribed  for  numerical 
reasons  only  in  order  to  remove  the  inevitable  accu- 
mulation of  enstrophy  on  the  grid  scale.  This  is  the 
reason  why  higher  order  operators  such  as  the  bihar- 
monic form  of  friction  have  traditionally  been  favored 
in  eddy-resolving  or  eddy-permitting  numerical  simula- 
tions ( Holland , 1978;  Bryan  and  Holland , 1989;  Smith 
et  al.,  2000).  Higher  order  operators  remove  numerical 
noise  on  the  grid  scale  and  leave  the  larger  scales  mostly 
untouched  by  allowing  dynamics  at  the  resolved  scales 
of  motion  to  dominate  the  subgrid-scale  parameteriza- 
tion ( Griffies  and  Hallberg,  2000). 

In  addition  to  numerical  closure,  the  viscosity  oper- 
ator can  also  be  a parameterization  of  smaller  scales. 
One  of  the  most  difficult  tasks  in  defining  the  param- 
eterization is  the  specification  of  the  Reynolds  stresses 
in  terms  of  only  the  resolved  scales’  velocities  [see  Ped- 
losky  (1979)  for  a review]  and  the  common  practice  has 
been  to  assume  that  the  turbulent  motion  acts  on  the 
large  scale  flow  in  a similar  manner  as  molecular  viscos- 
ity. However,  the  resulting  Laplacian  form  of  dissipa- 
tion removes  both  kinetic  energy  and  enstrophy  over  a 
broad  range  of  spatial  scales,  and  its  use  in  numerical 
models  in  general  implies  less  energetic  flow  fields  than 
in  cases  with  more  highly  scale-selective  dissipation  op- 
erators. In  order  to  assess  the  impact  of  the  dissipation 
operator  on  the  Gulf  Stream  system,  several  sensitivity 
experiments  were  performed  with  MICOM  using  Lapla- 
cian and  biharmonic  operators  for  the  viscosity  in  the 
momentum  equations. 

The  mean  SSH  of  two  simulations  performed  with 
two  different  magnitudes  of  the  biharmonic  viscosity 
coefficient  are  displayed  in  Figs.  3 and  4,  respectively 
(COADS-forced  run).  When  a relatively  small  value  of 
the  biharmonic  viscosity  coefficient  is  used  (see  caption 
of  Fig.  3 for  details),  the  western  boundary  current  is 
seen  to  separate  from  the  coast  early  at  the  Charleston 
bump  before  Cape  Hatteras  (Fig.  3).  A similar  result 
was  observed  with  the  1/10°  Los  Alamos  Parallel  Ocean 
Model  (POP)  during  the  spin-up  phase  in  which  both 


VISCOSITY  PARAMETERIZATION  AND  THE  GULF  STREAM  SEPARATION 


41 


Figure  3.  1-year-mean  model  SSH  field  with  a bihar- 
monic viscosity  operator;  A4  = max  [.1  Arc4  x deforma- 
tion tensor,  VbAx3],  with  Vd  = 1 cm/s. 


Figure  4.  1-year-mean  model  SSH  field  with  a bihar- 
monic viscosity  operator;  A4  = max  [Ax4  x deformation 
tensor,  VbAx3],  with  Vp  = 2 cm/s. 

the  viscosity  and  diffusion  had  to  be  increased  by  a 
factor  of  3 in  order  to  eliminate  this  feature  ( Smith 
et  al. , 2000).  An  increase  in  the  magnitude  of  the  bi- 
harmonic  viscosity  operator  in  MICOM  did  indeed  also 
eliminate  the  early  detachment  seen  in  Fig.  3,  but  it  also 
led  to  the  establishment  of  a permanent  eddy  north  of 
Cape  Hatteras  (Fig.  4).  This  eddy  results  from  a se- 
ries of  warm  core  (anticyclonic)  rings  that  propagate 
westward,  collide  with  the  western  boundary,  and  are 
only  weakly  dissipated  by  the  biharmonic  viscosity  op- 
erator. This  behavior  is  reminiscent  of  other  simula- 
tions performed  with  biharmonic  dissipation  ( Smith  et 


al. , 2000).  The  fact  that  this  permanent  eddy  only  ap- 
pears with  biharmonic  operators  seems  to  indicate  an 
incorrect  representation  of  the  eddy/mean  flow  and/or 
of  the  eddy/topography  interactions,  possibly  because 
of  the  scale  selectiveness  of  the  higher  order  operator 
that  allows  features  that  are  marginally  resolved  by  the 
grid  spacing.  In  all  simulations,  the  grid  spacing  is  such 
that  both  the  inertial  and  the  viscous  boundary  layers 
are  resolved  (very  well  for  the  inertial  and  minimally 
for  the  viscous). 


Figure  5.  1-year-mean  model  SSH  field  with  a Lapla- 
cian  viscosity  operator;  A2  = max  [.1  Ax2x  deforma- 
tion tensor,  VdAx],  with  Vp  — 1 cm/s. 

The  mean  SSH  of  the  simulation  performed  with 
the  Laplacian  viscosity  operator  is  displayed  in  Fig.  5 
(CO ADS- forced  run).  The  magnitude  of  the  Laplacian 
viscosity  coefficient  is  the  minimum  value  needed  for  nu- 
merical stability.  In  that  simulation  (Fig.  5),  the  Gulf 
Stream  separates  well  from  the  coast,  but  does  not  pen- 
etrate further  than  the  New  England  Seamounts. 

Overall,  neither  the  Laplacian  nor  the  biharmonic 
viscosity  operator  alone  provide  satisfactory  results  re- 
garding the  Gulf  Stream  system  behavior.  With  the 
biharmonic  operator,  eddies  are  found  to  retain  their 
structure  for  longer  periods  of  time  than  with  a Lapla- 
cian operator,  but  with  undesirable  effects  on  several 
features  of  the  large  scale  circulation.  With  the  Lapla- 
cian operator,  the  western  boundary  current  and  its 
separation  are  well  represented,  but  with  a weak  pene- 
tration of  the  Gulf  Stream. 

With  a Laplacian  (harmonic)  dissipation  operator, 
the  evolution  of  a wave  c(t)etkx  is  damped  exponentially 
with  a spin-down  time  t2  = AJ 1 (^sin  (^f5))  . In 

the  case  of  a biharmonic  operator,  the  spin-down  time  is 
T4  = Al 1 (^sin  . For  comparison  purposes, 
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Figure  6.  Laplacian  and  biharmonic  decay  time  scale 
as  a function  of  the  wavelength  k for  values  of  the  dif- 
fusive velocity  Vd  — -5  and  1 cm/s,  respectively. 


constant  harmonic  and  biharmonic  viscosity  coefficients 
can  be  expressed  as  a function  of  a diffusive  velocity 
Vd  and  the  grid  spacing  Ax  as  A2  = Vd  Ax  and  A4  = 
VdAx3,  respectively.  Examples  of  spin-down  times  for 
both  operators  are  given  in  Fig.  6 for  the  average  grid 
spacing  of  the  MICOM  simulations  (6  km).  For  the 
same  diffusive  velocity,  the  biharmonic  operator  more 
strongly  selects  the  small  scales  to  dissipate  and  leaves 
the  large  scales  relatively  untouched. 

The  Laplacian  experiment  of  Fig.  5,  when  contrasted 
to  the  biharmonic  experiments  of  Figs.  3 and  4,  sug- 
gests that  some  damping  of  the  larger  scales  is  neces- 
sary for  a reasonable  western  boundary  current  behav- 
ior. The  best  separation/penetration  results  were  ob- 
tained in  the  COADS-forced  and  the  ECMWF-forced 
runs  shown  in  Figs.  1 and  2 in  which  the  viscosity  op- 
erator was  prescribed  as  a combination  of  the  bihar- 
monic and  Laplacian  operators.  The  main  motivation 
for  combining  the  two  operators  (see  caption  of  Fig.  1 
for  details)  was  to  be  able  to  retain  the  scale  selective- 
ness of  the  biharmonic  operator  and  to  provide  some 
damping  at  the  larger  scales  [performed  in  this  case  by 
the  Laplacian  operator  for  k greater  than  80  km  (Fig. 
6)].  This  allowed  us  to  reduce  the  magnitude  of  the 
Laplacian  coefficent  A%  by  50%  and,  at  the  same  time, 
ensure  numerical  stability  with  an  effective  damping  of 
the  smaller  scales  via  the  biharmonic  operator  (Fig.  6). 
When  combined,  the  individual  diffusive  velocity  Vd 
specified  for  each  operator  is  smaller  than  the  mimi- 
mum  value  that  is  needed  for  numerical  stability  when 
only  one  of  the  operators  is  specified. 


4.  Summary  and  discussion 

These  results  appear  to  suggest  that,  in  a realistic 
setting,  even  with  such  a fine  grid  spacing,  the  mod- 
eled large  scale  ocean  circulation  is  strongly  dependent 
upon  the  choices  made  for  the  viscosity  operators.  Fur- 
thermore, it  appears  that  the  cascade  of  energy  from 
the  small  scales  to  the  larger  scales  may  not  take  place 
as  anticipated  and  that  some  large  scale  information 
is  needed  for  a proper  representation  of  the  western 
boundary  current.  In  the  experiments  described  in  this 
paper,  the  latter  is  taking  place  via  the  Laplacian  vis- 
cosity operator.  Hyperviscosity  (V2n  operator  with 
n > 2)  is  often  used  in  numerical  simulations  of  tur- 
bulent flows  to  extend  the  range  of  the  inviscid  iner- 
tial cascade.  It  has,  however,  been  argued  that  it  may 
also  contribute  non-trivial  spurious  dynamics  ( Jimenez , 
1994).  While  it  can  be  firmly  stated  that  a resolution  of 
1/10°  is  sufficient  for  the  Gulf  Stream  to  separate  from 
the  coast  at  Cape  Hatteras  ( Paiva  et  al.,  1999;  Hurlburt 
and  Hogan,  2000;  Smith  et  al.,  2000),  it  is  not  yet  clear 
what  is  the  optimal  resolution  for  a correct  Gulf  Stream 
penetration  and  variability.  A four-fold  increase  in  reso- 
lution from  1/16°  to  1/64°  with  the  Laplacian  operator 
in  the  hydrodynamic  (i.e.  no  thermal  forcing)  Navy 
Layered  Ocean  Model  (NLOM  - Hurlburt  and  Hogan, 
2000)  brought  the  SSH  variability  to  observed  levels 
without  altering  the  pattern  of  the  large  scale  circu- 
lation. While  numerical  simulations  at  the  above-noted 
resolutions  are  becoming  more  common,  they  still  de- 
mand the  latest  in  computing  facilities.  A four-fold  in- 
crease in  resolution  for  the  thermodynamically  forced 
models  cannot  be  realistically  implemented  with  the 
present  computer  resources.  Thus,  further  evaluation 
of  the  impact  of  various  dissipation  operators  on  the 
large  scale  circulation  should  be  pursued. 
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